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Simulations of binary black holes systems using the Spectral Einstein Code (SpEC) are done on a 
computational domain that excises the regions surrounding the black holes. It is imperative that the 
excision boundaries are outflow boundaries with respect to the hyperbolic evolution equations used 
in the simulation. We employ a time-dependent mapping between the fixed computational frame 
and the inertial frame through which the black holes move. The time-dependent parameters of the 
mapping are adjusted throughout the simulation by a feedback control system in order to follow 
the motion of the black holes, to adjust the shape and size of the excision surfaces so that they 
remain outflow boundaries, and to prevent large distortions of the grid. We describe in detail the 
mappings and control systems that we use. We show how these techniques have been essential in the 
evolution of binary black hole systems with extreme configurations, such as large spin magnitudes 
and high mass ratios, especially during the merger, when apparent horizons are highly distorted and 
the computational domain becomes compressed. The techniques introduced here may be useful in 
other applications of partial differential equations that involve time-dependent mappings. 

PACS numbers: 04.25.D-, 04.25.dg, 02.70.Hm 



I. INTRODUCTION 

Feedback control systems are ubiquitous in technolog- 
ical applications. They are found, for example, in ther- 
mostats, autopilots, chemical plants, and cruise control 
in automobiles. The purpose of a control system is to 
keep some measured output (such as the temperature in 
a room) at some desired value by adjusting some input 
(such as the power to a furnace). 

In the last few years, feedback control systems have 
also found applications in the field of numerical relativity, 
particularly in simulations of black hole binary systems 
that employ spectral methods and excision techniques [ - 
5]. 

Black hole excision is a means of avoiding the physi- 
cal singularities that lurk inside black holes. The idea is 
to solve Einstein's equations only in the region outside 
apparent horizons, cutting out the region inside the hori- 
zons. The boundaries of the excised regions are called 
excision boundaries. Causality ensures that the excision 
boundaries and the excised interiors cannot affect the 
physics of the exterior solution, and an appropriate hy- 
perbolic formulation of Einstein's equations [G, 7] can en- 
sure that gauge and constraint-violating degrees of free- 
dom also do not propagate out of the excised region. 

Excision is straightforward for black holes that remain 
stationary in the coordinates that are used in the simu- 
lation, but excision becomes more complicated when the 
black holes move or change shape. For numerical meth- 
ods based on finite differencing, the excision boundaries 
can be changed at every time step by activating or de- 
activating appropriate grid points and adjusting differ- 
encing stencils [8-15] . However, for spectral numerical 
methods, there is no equivalent of deactivating individ- 
ual grid points; instead, spectral methods are defined in 



finite extended spatial regions with smooth boundaries. 
The nearest equivalent to the finite-difference excision 
approach would be interpolating all variables to a new 
slightly offset grid at every time step, which would be 
computationally expensive. Therefore, spectral numer- 
ical methods need a different approach to reconcile the 
need for moving black holes with the need for a fixed 
excision boundary inside of each black hole. 

The solution [1] to this problem adopted by our group 
makes use of multiple coordinate systems. We call "iner- 
tial coordinates" those coordinates that asymptotically 
correspond to an observer at rest; in these coordinates 
the black holes orbit each other, have distorted shapes, 
and approach each other as energy is lost to gravitational 
radiation. Spectral methods are applied in another coor- 
dinate system, "grid coordinates" , in which the excision 
boundaries are spherical and stationary. We connect grid 
coordinates with inertial coordinates by means of an an- 
alytic mapping function M that depends on some set of 
time-dependent parameters X(t). These parameters must 
be continually adjusted so that each spherical, stationary 
grid-frame excision boundary is mapped to a surface in 
the inertial frame that follows the motion and the shape 
of the corresponding black hole as the system evolves. It 
is this adjustment of each parameter A(t) that is accom- 
plished by means of feedback control systems, one control 
system per parameter. 

In this paper, we describe in detail the mapping func- 
tions and the corresponding feedback control systems 
that we use to handle black hole excision with spec- 
tral methods. Some earlier implementations have been 
described previously [1, 2, 4, 16, 17] , but there have 
been many improvements that now allow spectral exci- 
sion methods to produce robust simulations of binary 
black hole systems, including those with unequal masses, 
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high spins, and precession. We describe these improve- 
ments here. In Sec. II we review control theory and 
present simple examples of control systems. In Sec. Ill 
we discuss the implementation of control systems in the 
SpEC [18] code. Section IV details the coordinate map- 
pings used in SpEC and the feedback control systems 
used to control them. In Sec. VI we describe applica- 
tions of control systems in SpEC besides the ones used 
to adjust map parameters. We summarize in Sec. VII. 



II. CONTROL THEORY 

To motivate control theory, we begin with a simple 
example: cruise control in an automobile. Suppose we 
wish to control the speed v of a car that is driving up an 
incline of small angle 9. The equation of motion for this 
system is 



dv 
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m m 



(1) 



where m is the mass of the car, g is the gravitational 
acceleration, r\ is a drag coefficient, and F is a force sup- 
plied by the car's engine. We wish to determine F so as 
to cause the car to maintain a speed of v = vq, even if 
the angle 9 changes as the car climbs the incline. To do 
this, we choose 



F 



= K P Q + K T / Q{t) dr, 



(2) 



where Kj and Kp are constants, and where Q = Vq — v 
is the quantity that we wish to drive to zero. We call Q 
the control error. Differentiating Eq. (1) and substituting 
Eq. (2) yields 
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dt 2 
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which is the equation for a damped, forced harmonic os- 
cillator. By choosing Kp to produce critical damping, 
and choosing K\ to set a timescale, this choice will drive 
v towards vq as desired. 

The basic structure of a control system is easily un- 
derstood as a feedback loop (see Fig. 1). The simulation 
(e.g., a binary black hole simulation, or a car driving up 
an incline) produces some measure of error, Q(t), that de- 
fines the deviation from some desired target value. This 
error acts as the input for the control system, which then 
computes a control signal U(t) (e.g., the derivative of a 
map parameter, or the force supplied by a car's engine) 
that will minimize the error Q(t) when fed back into the 
simulation. 

A simple and effective way to compute U(t) is to make 
it a linear combination of the error, Q(t), and integrals 
and/or derivatives of the error. The term proportional 
to the error acts to reduce the deviation from the de- 
sired value, the terms proportional to derivatives of the 
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FIG. 1: A generic control circuit. The simulation outputs a 
measure of error, Q, which is used by the control system. The 
control system then outputs a signal, U, which changes the 
behavior of the simulation. 



error act to oppose rapid deviations, and the terms pro- 
portional to the integrals act to reduce any persistent 
deviation or offset that accumulates over time. In the 
cruise control example, we used a proportional and an 
integral term only in Eq. (2). 

We now turn to another example of a control system 
that is more closely related to the way we use control 
systems in binary black hole simulations. Consider two 
coordinate systems, (x,t) and (x,t), that are related by 
the map 



x 
t 



x 
i. 



V(t)t, 



(4) 
(5) 



We wish to control the parameter V(t) so that a wave 
f{x, t) = f(x — ct) that propagates at speed v = c in 
the (x,t) coordinates will propagate at some arbitrary 
desired speed v d in the (x, t) coordinates. According to 
Eq. (4), the speed of the wave in the (x,t) coordinates is 
v = c— V(t). Therefore we define a control error Q to be 



Q = c- V(t) - v dl 



(6) 



and we construct a control system that drives this control 
error to zero. If we choose the control signal U(t) to be 
d 2 V/dt 2 , then the simplest feedback loop that can be 
constructed uses only a term proportional to the error 
and amplified by a "gain" Kp. Then the evolution of 
V(t) is given by 



d 2 V 
dt 2 



K P Q = K P [c- V(t) - v d ] . 



The solution to this equation is of the form 

V(t) — c — Vd + Ax sin(crf) + A 2 cos(crf), 



(7) 



(8) 



where a :— ^/Kp, which is oscillatory for Kp > and 
divergent for Kp < 0. Thus we see that adding only 
a proportional term to the control signal U(t) is insuffi- 
cient, since it does not reduce the control error Q. 

However, if we add a derivative term to the feedback 
equation, 



K P Q + K D 



d 2 V 
dt 2 

then the solution is of the form 



dQ 

dt ' 



(9) 



V(t) 



c—v d +cxp 



i KDt [Bx sin(/3<) + B 2 cos{pt)] , (10) 
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where ft := ^^f&Kp — K^. This solution is stable with 
an exponentially damped envelope when 4Kp > K^, 
which will cause v — ¥ Vd as t —¥ oo. 

Notice that this control system allows us to choose a 
V(t) such that v is the opposite sign of c and the wave is 
left-going in the (x, t) frame instead of right-going. This 
behavior can be seen in the example in Fig. 2. 
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FIG. 2: The speed in (x, t) coordinates, v, is plotted for a 
family of gains Kp with Kp> — 1 fixed. The wave velocity is 
c = —0.2 and the desired velocity is Vd = 0.5 (dashed line). 
The controller turns on at t = 2. For gains in the stable 
region, where Kp > 0.25, v settles down to Vd- One over- 
damped situation with Kp = 0.01 is plotted for comparison 
(dotted line). 

The overdamped solution (dotted line in Fig. 2) has 
a persistent offset that can be ameliorated by adding an 
integral term to the feedback equation, as was done in the 
cruise control example. In principle we could continue to 
add more terms, but in practice, it is usually sufficient 
to use a PID (proportional-integral-derivative) controller, 
which has terms proportional to the control error, its 
integral, and its derivative. When the underlying system 
is unknown this is the best controller to use [19]. 



III. CONTROL SYSTEMS IN SPEC 

Control systems are used in SpEC for several purposes. 
The most important is their role in handling moving, ex- 
cised black holes in a spectral evolution method. We 
use a dual-frame method [1] in which the grid is fixed in 
some coordinates {t,x l ) but the components of dynami- 
cal fields are expressed in a different coordinate system 
it, x l ). We call (t, x l ) the grid coordinates and it, x l ) 
the inertial coordinates. Figure 3 shows an example do- 
main decomposition in both coordinate systems. The 
two coordinate systems are connected by a map M that 
depends on time-dependent parameters X(t). The exci- 
sion boundaries are exactly spherical and stationary in 
grid coordinates. In inertial coordinates, the apparent 
horizons move and distort as determined by the solution 
of Einstein's equations supplemented by our gauge con- 



ditions. The parameters A(i) need to be controlled so 
that the excision boundaries in the inertial frame follow 
the motion and shapes of the apparent horizons. We use 
control systems to accomplish this. 




FIG. 3: a) Computational domain in grid coordinates; the 
black hole centers are at rest and the excision boundaries 
are spherical, b) Same domain in inertial coordinates near 
merger; the excision boundaries move and distort to track 
the apparent horizons. 

The particular maps that we use will be discussed in 
Sec. IV. In this section we will describe how we construct 
the control system for a general parameter A(t), includ- 
ing how we define the relationship between the control 
error Q(t) and the control signal U(t), how we smooth 
out noise in the control system, and how we dynamically 
adjust the feedback parameters and timescales. 

A. Definition of control errors and control signals 

We represent a general time-dependent map parameter 
A(t) as a polynomial in time with a piecewise constant 
Nth derivative: 

N 1 

A(t) = y)- f (*-t i ) n A? J for ti<t<t i+1 , (11) 

n=0 

where for each time interval ti < t < ij+i the quantities 
A" are constants. 

At the beginning of each new time interval ti, we 
set the constants A™ in Eq. (11) as follows. First for 
n = N we set A^ = Ufa), where U{t) is the con- 
trol signal defined in detail below. For n < N we set 
A" = d n X/dt n \ t=t , where the derivative is evaluated at 
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the end of the previous time interval. In this way all 
the derivatives of X(t) except the Nth derivative are con- 
tinuous across intervals. The goal will be to compute 
the control signal U(t) so as to drive the map parameter 
X(t) to some desired behavior. Before we describe how 
to compute the control signal U(t), we first discuss the 
control error Q, which will be used in the computation 
of [/(f). 

To appropriately define the control error Q, one must 
answer the question of how a small change in the map 
parameter corresponds to a change in the observed vari- 
ables. If the control error is defined to be too large, then 
the controller will consistently overshoot its target, po- 
tentially leading to unstable behavior; conversely, if the 
control error is defined to be too small, then the controller 
may never be able to reach its target value. 

If there exists a target value of X(t), call it Atarget, that 
does not depend on the map but may depend on other 
observable quantities in the system (call them A, B,. . . ), 
then wc define 



Q = Atarget {A,B,...) - A. 



(12) 



The goal is to drive Q to zero and thereby drive A to 

Atarget • 

If instead, as often happens for nonlinear systems, the 
target value of A depends on A itself, even indirectly, then 
we define Q differently using a generalization of Eq. (12): 
we require that A attains its desired value when Q — > 0, 
and we require that 



dQ 
dX 



1 + 0(Q). 



(13) 



Note that in either Eq. (12) or Eq. (13), Q could in prin- 
ciple be multiplied by an arbitrary factor; however, if 
this were done, that factor would need to be taken into 
account in the computation of the control signal U(t) be- 
low. Without loss of generality we assume no additional 
scaling. 

In the case of several map parameters X a (t) with cor- 
responding Q a , where a is an index that labels the map 
parameters, the desired value of some A a may depend on 
a different map parameter At,. In this case, we generalize 
Eq. (13) and define the Q a associated with A Q as follows: 



dQ a 
dX b 



0(Q), 



(14) 



where 5 a b is a Kronecker delta. This criterion ensures 
that we can treat each A a independently when all control 
errors are small. A way of understanding Eq. (14) is to 
consider a set of Q' a that are obtained via Eq. (13) with- 
out regard to coupling between different A Q . Then a set 
of Q a satisfying Eq. (14) can be obtained by diagonaliz- 
ing the matrix dQ' a /dX D . In the remainder of this section 
we assume that if there are multiple map parameters A Q , 
the corresponding control errors Q a satisfy Eq. (14). We 
therefore drop the a indices and write equations for U(t) 



in terms of a single Q(t) satisfying Eq. (13) that repre- 
sents the control error for a single X(t). 

The control error Q(t) is a function of several observ- 
ables. In the case of the mapping functions that are 
designed to move and distort the excision boundaries to 
follow the motion and shapes of the apparent horizons, 
Q(t) is some function of the current position or shape of 
one or more apparent horizons. The precise definition of 
Q(t) is different for each map parameter, and depends 
on the details of how each map parameter couples to the 
observables. We will discuss the control error Q for each 
of the map parameters in Sec. IV. But it is not necessary 
to know the exact form of the control error in order to 
compute the control signal U(t); it suffices to know only 
that U (ti) is equal to Af in Eq. (11), and that the control 
error Q obeys Eq. (13). 

We now turn to the computation of the control signal 
U(t). There is some flexibility in the control law deter- 
mining U(t), so long as key feedback mechanisms are in 
place (as shown in Sec. II). In SpEC, we use either a 
standard PID controller, 



U(t) =a Q Q(t) dt + at Q(t) + a 2 



dQ 

dt ' 



or a special PD controller, 



K 



k=0 



(15) 



(16) 



where typically K = 2. 

We set the constants so that the system damps Q 
to zero on some timescale Td that we choose. We assume 
that Td is longer than the interval tj+i — f defined in 
Eq. (11), so that we can approximate Q(t) and X(t) as 
smooth functions rather than as functions with piecewise 
constant iVth derivatives. Under this assumption, we 
write 



U(t) 



d N X/dt N , 
-d N Q/dt I 



(17) 



where in the last line we have used Eq. (13), we have as- 
sumed that Q is small, and we have neglected the time de- 
pendence of other parameters besides A that enter into Q 
under the assumption that the control system timescale 
is shorter than the timescales of the quantities that we 
want to control. For the PID controller and N = 2, com- 
bining Eqs. (15) and (17) yields 

" Q a I ' Q{t)dt + a l Q{t) + a 2 d Q. (18) 



dt 2 



If we choose a = 1/r?, a\ = 3/tJ, and a 2 = 3/t,j, then 
the solution to Eq. (18) will be exponentially damped on 
the timescale Td, 

-t/r d 



Q oc e 



(19) 



The same exponential damping holds for the PD con- 
troller, Eq. (16), for appropriate choices of the parame- 
ters cjfc. For K = 2 and N = 3, the parameters are 
identical to those in the PID case above. 
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B. Averaging out noise 

The PID controller, Eq. (15), is computed by measur- 
ing the control error Q, its time integral, and its time 
derivative. The PD controller, Eq. (16), may require 
multiple derivatives of Q{t) depending on the order K. 
Generally only Q, and not its derivatives or integrals, is 
available in the code at any given time step. The simplest 
way to compute the derivatives of Q is by finite differenc- 
ing in time, and the simplest way to compute the integral 
is by a numerical quadrature. 

We measure the control error Q at time intervals of 
length t to , where we choose r m < — £j. The mea- 
suring time interval r m is then used as the time step for 
hnite difference stencils and quadratures. 

The measured Q is typically a function of apparent 
horizon locations or shapes, and this measured Q includes 
noise caused by the hnite resolution of the evolution, and 
the finite residual and finite number of iterations of the 
apparent horizon finder. Taking a numerical derivative of 
Q amplifies the noise, and then this noise is transferred 
to the control signal via Eq. (15) or (16), and then to 
the map. If the noise amplitude is too large, the con- 
trol system will become unstable. The PID controller 
generally handles noise better than the PD controller for 
two reasons. First, each successive numerical derivative 
amplifies noise even further, so using only one derivative 
instead of two (or more) results in a more accurate con- 
trol signal. Second, the inclusion of an integral term acts 
to further smooth the control signal. 

In some cases, however, the noise in Q can be prob- 
lematic even for the PID controller. In these cases we 
implement direct averaging of the control error in one of 
two ways: 1) we perform a polynomial fit of order N to 
the previous M measurements of the control error, where 
M > N, or 2) we perform an exponentially- weighted av- 
erage, with timescale T avg , of all previous control error 
measurements and their derivatives and integrals. 

The latter method is implemented by backwards-Euler 
differencing of the following ODE: 

^=F(i)-F avg /r avg , (20) 

where F(t) is the measured value of Q, d n Q/dt n , or 
J Q dt at a particular time, and F avg is the averaged value 
of this quantity, i.e., the value that is substituted into 
Eq. (15) or (16). The backwards-Euler differencing takes 
the form 

pk+i _ pfe 

aVS x aVS = F(t k+1 ) - F^+Vravg. (21) 

> m 

We choose r avg to be less than the measurement timescale 
typically r avg ~ 0.25r m . 



C. Dynamic timescale adjustment 

In the previous section we have identified four 
timescales relevant for each control system. The first 
is the damping timescale Td] this describes how quickly 
the control error Q(t) falls to zero, and therefore how 
quickly the map parameter X(t) approaches its desired 
value. The second timescale is the control interval Ati := 
ti+i—ti, which represents how often the Nth derivative of 
the function X(t) in Eq. (11) is updated. The third is the 
measurement timescale r m , which indicates how often the 
control error Q is measured. The fourth is the averaging 
timescale r avg , which is used to smooth the control error 
Q (and its derivatives and integrals) for use in comput- 
ing the control signal U(t). These timescales are not all 
independent; for example we have assumed Ati < T d hi 
deriving Eq. (18), and we have assumed T m < Ati so that 
we can obtain smooth measurements of the derivatives of 

Q. 

Because binary black hole evolutions are nonlinear dy- 
namic systems, we adjust the damping timescale, T d , 
throughout the simulation. We then set the three other 
timescales, Ati, T m, and r avg , based on the current value 
of the damping timescale Td, as we now describe. 

The timescale Td should be shorter than the timescale 
on which the physical system changes; otherwise, the con- 
trol system cannot adjust the map parameters quickly 
enough to respond to changes in the system. But if the 
timescale Td is too small, then the measurement timescale 
T m on which we compute the apparent horizon must also 
be small, meaning that frequent apparent horizon com- 
putations are needed; this is undesirable because comput- 
ing apparent horizons is computationally expensive. We 
would like to adjust Td in an automatic way so that it is 
relatively large during the BBH inspiral, decreases during 
the plunge and merger, and increases again as the rem- 
nant black hole rings down. In the canonical language of 
control theory, we would like to implement "gain schedul- 
ing" [20], tuning the behavior of the control system for 
different operating regimes. 

We do this as follows: For all map parameters (except 
the £ — component of the horizon shape map, which is 
treated differently; see Sec. VD), the damping timescale 
is a generic function of Q and Q, i.e., the error in its 
associated map parameter and its derivative. Whenever 
we adjust the control signal U(t) at interval ti, we also 
tune the timescale in the following way: 

Td +1 = (22) 

where typically 

f 0.99, if Q/Q > -1/2t 4 and \Q\ or \Q Td \ > Qf ax 
P = < 1.01, if \Q\ < Qf™ and |Qr d | < Qf in 
[ 1, otherwise. 

(23) 

Here Q^ ln and Q^ ax are thresholds for the control er- 
ror Q. The idea is to keep Q < Q^ ax so that the map 
parameters are close to their desired values, but to also 



() 



keep Q > Q^ 1 " 1 because an unnecessarily small Q means 
unnecessarily small timescales and therefore a large com- 
putational expense (because the apparent horizon must 
be found frequently). For binary black hole simulations 
where the holes have masses Ma and Mb, we find that 
the following choices work well: 



}Max 



2 x 10" 



M A /M B + M B /M A 



Q 



Min 



(Max 



(24) 
(25) 



Once we have adjusted the timescale for every con- 
trol system, we then use these timescales to choose the 
times ti+i in Eq. (11) at which we update the polynomial 
coefficients A™ in the map parameter expression A(t), 



Mi := U +1 -U = a d min(T d ), 



(26) 



where typically ay — 0.3, and the minimum is taken over 
all map parameters (except for the £ = component of 
the horizon shape map). This ensures that the coeffi- 
cients A™, are updated faster than the physical system is 
changing, and faster than the control system is damping. 
For ad too large, we find that the control system becomes 
unstable. 

We also use the timescale to choose the interval r TO 
at which we measure the control error. For many map 
parameters, the associated control error is a function of 
apparent horizon quantities, which is why we desire r m 
to be as large as possible. But a certain number of mea- 
surements are needed for each Ati so that the control 
signals, defined in Eqs. (15) and (16), are sufficiently ac- 
curate and the control system is stable. We choose 



(27) 



where a m is typically between 0.25 and 0.3. In other 
words, we measure the control error three or four times 
before we update the control signal. This also sets the 
averaging timescale (see Sec. Ill B) as we typically choose 

7"ave: 0.25T m . 



IV. CONTROL SYSTEMS FOR MAPS 

In a SpEC evolution, we transform the grid coordi- 
nates, x % , into inertial coordinates, x 1 , through a series 
of elementary maps as in [1-5]. Several maps have been 
added and many improvements have been made since 
their original introduction. The full transformation is 
x l = Mx 1 , where 



M, — M Translation ° -^Rotation ° -^Scaling 
°M S kew ° McntX ° Ms • 



I Shape ■ 



(28) 



Below we will describe each of these maps and how we 
measure the error in their parameters. Before we do this, 
however, we describe our domain decomposition and how 



we measure apparent horizons, because information from 
the grid and the horizons is used to determine the maps. 

In the grid coordinates x l , the domain decomposition 
looks like Fig. 4. There are two excision boundaries, A 
and B, which are spheres in grid coordinates. The grid- 
coordinate centers of these excision boundaries we will 
call C H , where H is either A or B. The excision bound- 
aries (and therefore C H ) remain fixed throughout the 
evolution. The purpose of many of the maps is to move 
the mapped centers C H := Ai (C H ) along with the centers 
of the apparent horizons. 

We measure the apparent horizons in an intermediate 
frame x l , which we call the distorted frame. This frame 
is connected to the grid frame by the map 



M 



Distorts 



M 



CutX 



oM 



Shape- 



(29) 



We will discuss exactly what M shape and .Mcutx are 
below, but for now we need only demand that M. Distortion 
has two properties: The first is that it leaves the centers 
of the excision boundaries invariant, i.e., 



C 



H 



M Distortion (C l H ) = C 



H- 



(30) 



The second property is that at each excision boundary H, 
-^Distortion leaves angles invariant. That is, if we define 
grid- frame polar coordinates (r# , Oh, </>#) centered about 
excised region H in the usual way, 



x a = + r H sin 9h cos (j>n, 
x 1 

„2 _ /~<2 



1 — C}f + r H sinflff sin^# , 



x z = Cff + r H cos6 H , 



(31) 
(32) 
(33) 



and similarly for polar coordinates (fn , &h, 4>h) centered 
about C\j in the x l frame, then the second property 
means that 



in 



Oh, 
<Ph, 



(34) 



on the excision boundary. Note that the index i on spatial 
coordinates ranges over (0, 1,2) in this paper. 

We find the apparent horizons in the frame x % by a 
fast flow method [21]. Each horizon is parameterized 
as a smooth surface in this frame, and is parameterized 
in terms of polar coordinates centered around C H : the 
radius of the horizon at each (9h, 4>h) is given by 



f H R (0H, <Ph) = > ' Sf m Y im (9 H , 4>h) 



H 



(35) 



where again the index H labels which excision boundary 
is enclosed by the apparent horizon. Here Yi m (0H,4>H) 
are spherical harmonics and S^ m are expansion coeffi- 
cients that describe the shape of the apparent horizon. 

We search for apparent horizons in the distorted frame 
x l rather than in the grid frame x l or the inertial frame 
x l because this simplifies the formulae for the control 
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systems. In particular, this choice decouples the con- 
trol errors of the shape map Qi m (defined in Sec. IV E, 
Eq. (76)) for different values of (£,m), and it decouples 
the errors Qe m from the control errors of the maps con- 
necting the distorted and inertial frames. 

For each surface H, we define the center of the appar- 
ent horizon 



(36) 



where the integrals are over the surface. In the code, it 
suffices to use the following approximation of Eq. (36), 
which becomes exact as the surface becomes spherical: 



& = C° H - v/sT^Re^) 
=C l H + v^Trlm^) 



(37) 
(38) 
(39) 



Here we have used the property C l H = C l H , Eq. (30). We 
distinguish the center C l H of the excision boundary from 
the center £]y of the corresponding apparent horizon. The 
former is fixed in time in the grid frame, but the latter 
will change as the metric quantities evolve. The purpose 
of several of the maps (namely M scaling, M Rotation, and 
■M Translation described below) is to ensure that — C l H 
is driven towards zero; i.e., to ensure that the centers of 
the excision boundaries track the centers of the apparent 
horizons. 

We now describe each of the maps comprising A4 and 
how we measure the error in their parameters. The order 
in which we discuss the maps is not the same as the order 
in which the maps are composed so that we can discuss 
the simplest maps first. 



A. Scaling 

The scaling map, M. Scaling, causes the grid to shrink 
(and the excision boundaries to move closer together) in 
the inertial frame, thereby allowing the grid to follow the 
two black holes as their separation decreases because of 
gravitational radiation. 

This map transforms the radial coordinate with respect 
to the origin (i.e., the center of the outermost sphere in 
Fig. 4), such that the region near the black holes is scaled 
uniformly by a factor a and the outer boundary is scaled 
by a factor b, 



R h-> aR + (b - a)R 3 1 R 2 OB . 



(40) 



Here R is the radial coordinate, Rob is the radius at the 
outer boundary, a is a parameter that will be determined 
by a control system, and b is another parameter that will 
be determined empirically. 

In [ ], the scaling map was simply x l H ¥ a{t)x l , which 
is recovered by Eq. (40) if b = a. In that case, as the 



black holes inspiral together and a decreases, the outer 
boundary of the grid decreases as well. For long evo- 
lutions, the outer boundary decreases so much that we 
can no longer extract gravitational radiation far from the 
hole. The addition of b to the map alleviates this diffi- 
culty, allowing the motion of the outer boundary to be 
decoupled from the motion of the holes. We choose b by 
an explicit functional form 



b(t) = 1 - 10- b *7(2500 + r) 



(41) 



which is designed to keep the outer boundary from 
shrinking rapidly, but to allow the boundary to move 
inward at a small speed, so that zero-speed modes are 
advected off the grid (and thus need no boundary condi- 
tion imposed on them). We find that this reduces errors 
at the outer boundary. 

We control the scale factor a of this map using the 
error function 



where 



Q a = a{dx"-1), 



dx< = €a Cb 



(42) 



(43) 



We assume that the separation vector C\ — C l B in the 
grid frame is parallel to the ic-axis. The idea is that the 
distorted-frame separation of the horizon centers along 
the x-axis is driven to be the same as the separation of 
the excision boundary centers. 

To show that the control system for a obeys Eq. (13), 
we consider the change in Q a under variations of a with 
all other maps held fixed, and with the inertial-coordinate 
centers of the horizons held fixed. 1 This means that 
the distorted-frame centers of the horizons £ l H appear- 
ing in Eq. (43) change with variations of a. The second 
term in Eq. (40) is small when evaluated near the horizon 
where (R/Rob) 2 <C 1, so we can write the action of the 
scaling map as 



and therefore under variations 5a, 

= Se H =C H Sa + a6i i H . 



(44) 



(45) 



Taking variations of Eq. (42), using Eq. (45) to substitute 
for S^ l Hl and noting that the excision-boundary centers 
C 1 ^ are constants and do not vary with a, we obtain 



SQ a = —5a, 



(46) 



We consider inertial-framc quantities to be fundamental and 
determined by the solution of Einstein's equations plus gauge 
conditions, and therefore independent (modulo numerical trun- 
cation error) of the numerical grid and of the maps that we use 
to construct, move, and distort that grid. 
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so that Eq. (13) is satisfied. 

To verify that the more general decoupling equation, 
Eq. (14), is satisfied for the scaling map, we must show, 
in addition to Eq. (46), that the quantity dx° defined in 
Eq. (43) is invariant under the action of the other maps, 
at least in the limit that all the control errors Q are 
small. We consider each map in turn. The translation 
map moves both apparent horizons together, so it leaves 
da; invariant. Changes in the rotation map parameters 
will change dx° only by an amount proportional to the 
control errors of the rotation map (see Eqs. (52) and (53) 
below). The skew map (below) leaves the centers of the 
excision boundaries invariant. Because the intent of the 
rotation, translation, and scaling maps is to drive the 
centers of the apparent horizons toward the centers of 
the excision boundaries, this means that the skew map 
changes dx° by an amount proportional the control errors 
of the rotation, translation, and scaling maps. Finally, 
the shape and CutX maps connect the distorted and grid 
frames, so they cannot affect dx°. 



B. Rotation 

The rotation map, A^Rotation, is a rigid 3d rotation 
about the origin that tracks the orbital phase and pre- 
cession of the system, 



cos ■& cos tp — sin $ cos sin tp 
— I sin z9 cosy cost? sin?? sin y | .<•• 
— sin tp cos y 



(47) 



The pitch and yaw map parameters (tp, $) are controlled 
so as to align the line segment connecting the apparent 
horizon centers with the distorted-frame rc-axis. Note 
that the map parameters (tp, $) are functions of time, 
and are not to be confused with the polar coordinates 
(4>Ht9h) centered about each excision boundary. Then 
the control error is given by 



e B 



e A -e B 
a - & 

(£°-|°) cosy' 



(48) 
(49) 



In the case of extreme precession where <p — > tt/2, these 
equations are insufficient because Q$ diverges. Our so- 
lution is to use quaternions, which avoid this singularity 
(for a complete discussion, see Ossokine 2012, in prepa- 
ration). 

The control errors Q$ and Q v obey Eq. (14). To show 
this, consider variations of $ and tp with other maps held 
fixed, and with the inertial-coordinate centers £jy of the 
horizons held fixed. The rotation map, Eq. (47), implies 
that under these variations, 



where 1Z is the matrix in Eq. (47) . Multiplying this equa- 
tion by 7?. -1 we obtain 



5e H = -n- l (8ii)C H , 



which yields 

dtp 
3d 




cos tp 
cosy — sin tp | 
sin tp 



(51) 

(52) 
(53) 



We can now verify Eq. (14) for the special case where the 
indices a and b in Eq. (14) are either d and tp. This result 
is obtained in a straightforward way by differentiating 
Eq. (48) or (49) with respect to § or tp, and substituting 
Eqs. (52) or (53). 

In addition, the control errors Q v and Q$ are indepen- 
dent (to leading order in the control errors) of changes in 
the parameters of the other maps that make up M. : Vari- 
ations of Eqs. (48) and (49) with respect to the scaling 
map parameter a are zero, because both the numerator 
and denominator of Q v and Q-o scale in the same way 
with a. Similarly, Q v and Qs are independent of the 
translation map, since both apparent horizons are trans- 
lated by the same amount. The skew map can change 
Qtp and Q§, but only by an amount proportional to con- 
trol errors, because the skew map leaves C l H invariant 
and other maps ensure that £ l H are close (within a con- 
trol error) to Cjj. The shape and CutX maps cannot 
affect Q v and Q$ because those maps connect the grid 
and distorted frames (and therefore they cannot change 
the distorted- frame horizon centers £]y). 



C. Translation 

The translation map, M Translation, transforms the 
Cartesian coordinates, x l , according to 



x i y x 



f(R)T l 



(54) 



(50) 



where f(R) is a Gaussian centered on the origin with 
a width set to the outer boundary radius and T l (t) are 
translation parameters that are adjusted by a control sys- 
tem. 

The translation map moves the grid to account for any 
drift of the "center of mass" of the system (as computed 
assuming point masses at the apparent horizon centers) 
in the inertial frame. This drift can be caused by mo- 
mentum exchange between the black holes and the sur- 
rounding gravitational field [22, 23], or by anisotropic ra- 
diation of linear momentum to infinity. This is the third 
map (the other two are rotation and scaling) that drives 
the centers of the apparent horizons toward the centers 
of the excision boundaries. The apparent horizon cen- 
ters t; l A and !; l B represent six degrees of freedom: one is 
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fixed by the scaling map, two by rotation, and three by 
translation. 

The control errors will be more complicated than for 
the other control systems because translation and rota- 
tion do not commute. We define the control errors Q l T 
for each of the translation directions i = 0, 1, 2 as 



Q'-i 



aK 



(55) 



where 1Z is the rotation matrix in Eq. (47), and V is some 
matrix yet to be determined, which may depend on the 
rotation parameters (<p>,$) and on the constants C B , but 
which may not depend on the translation parameters T l . 
This control error must have the property that Q l T — 
when £jj 
satisfies 



C l H , so our first restriction on V is that it 



= C B +V(C A -C B ). 



(56) 



To check Eq. (14), we note that near the black holes 
we can neglect the last term in Eq. (40) and we can use 
f(R) ~ 1 in Eq. (54), so that the apparent horizon cen- 
ters in the inertial and distorted frames are related by 

ijj = a- 1 ^ 1 ^ - a^n^T. (57) 

Inserting Eq. (57) into Eq. (55), we obtain 

Q^ = e B -T i -nvn- 1 (e A -e B )- (58) 

The second term in Eq. (58) is the only term that depends 
on the the translation parameter T l , so we have 



dQ l T 



(59) 



When varying map parameters, the inertial-frame hori- 
zon centers remain fixed, so the only other term in 
Eq. (58) that depends on map parameters is the last 
term, which depends on tp) because of the rotation 
matrices and because of the (f?, tp) dependence in V . We 
can therefore write the variation of Q l T with respect to 
($, tp) as 



dQ l T 
dW 



d 
dW 

d 
dW 



(lZVlZ~ 1 )alZ(C A - 



(60) 
(61) 



where W stands for either •& or tp, and where in the last 
line we have used Eq. (57). To obey Eq. (14), dQ^/dW 
must either be zero, or on the order of a control error Q. 
For i = 1,2, the quantity (£ A — £ B ) is proportional to the 
control error of the rotation map, Eqs. (48) and (49), so 
these terms can be neglected in Eq. (61). However, for 
* = 0, (C A — Cb) i s n °t proportional to a control error; 
instead (£ A — £,%) is driven to a constant finite value of 
C A — C B by the scaling control system. Therefore, in 
order to satisfy Eq. (14), we require 



d 
dW 



(KVK-^ll 



0. 



(62) 




FIG. 4: Two-dimensional projection of the domain decom- 
position near the two black holes A and B. Shown are bound- 
aries between subdomains. Each subdomain takes the shape 
of a spherical shell, a distorted cylindrical shell, or a distorted 
cylinder. Additional spherical-shell subdomains (not shown) 
surround the outer boundary of the figure and extend to large 
radius. This domain decomposition is explained in detail in 
the Appendix of [5] . The red, magenta, and blue surfaces are 
those for which the function /a {ta , Q a , 4>a) from Eq. (71) has 
a discontinuous gradient, as described in Sec. IV E. 



We find that we can satisfy both Eqs. (56) and (62) by 
choosing 



V = 



1 



c% 



-ch 

- C B tan if 
-C B tan tp 



C, 2 



C° B 



(63) 



Here we have again assumed that the separation between 
the centers of the excision boundaries is parallel to the 
C\ 



x-axis, i.e. 



and C\ 



(l 2 



D. Skew 

In Fig. 4, a prominent feature of the domain decom- 
position is a plane (a vertical line in the two-dimensional 
figure) that is perpendicular to the grid-frame x-axis and 
lies between the two excision boundaries A and B. We 
call this plane the "cutting plane" . 

The skew map, A^skow: acts on the distorted-frame 
coordinates, in which the cutting plane is still perpendic- 
ular to the x-axis. The skew map leaves the coordinates 
y, z unchanged, but changes the x-coordinate in order to 
give a skewed shape to the cutting plane, as shown in 
Fig. 5. Let x l c be the intersection point of the segment 



10 




FIG. 5: Example of the skew map in the mapped frame. 
Shown is a near-merger snapshot of an equal mass system 
where the holes have spin \ = 0.95 aligned with the orbital 
angular momentum. We see that an unskewed cutting plane 
(i.e., a vertical line) would be inadequate for this configura- 
tion because the cutting plane would intersect the excision 
boundaries. 



connecting the excision boundary centers and the cutting 
plane. The action of the skew map is denned as 

x°^x° + W tan [ 0J (*)] ' ( xJ - x c)> ( 64 ) 
i=i,2 

x 1 i ^ x j , for j = l,2 (65) 



FIG. 6: A diagram of the skew map in the mapped frame. The 
horizon surfaces are drawn in black and the excision surfaces 
are drawn in red. The dotted red line represents the segment 
connecting the excision surface centers, which is parallel to 
the a;-axis. The normal to either horizon at the intersection 
point with this segment is represented by the straight black 
lines. The cutting plane is represented by the (skewed) blue 
curve, and the normal to this plane at Xq is represented by the 
three green lines. There are three angles involved in the skew 
control system - the angle between the normal to the skewed 
cutting plane and the x-axis, Q J , and the angle between the 
normal to either horizon and the normal to the skewed cutting 
plane, O^. Skew control acts to minimize Q 3 H . We measure 
Q 3 H in the unmapped (distorted) frame, but is invariant 
under the skew map for W ~ 1. 



where W is a radial Gaussian function centered around 
ip, and the angles Q J (t) are the time-dependent map pa- 
rameters. For j = 1,2 the parameter J '(i) is the angle 
between the mapped and unmapped x-'-axis when pro- 
jected into the x 3 " 3 = x\j 3 plane. Note that the line 
intersecting x l c and parallel to the x-axis is left invariant 
by the skew map 2 . The width of the Gaussian W is set 
such that W is below machine precision at the innermost 
wave extraction sphere. This implies that the spherical- 
shell subdomains used to evolve the metric in the wave 
zone will not be affected by the skew map. 

The purpose of the skew map is to (as much as pos- 
sible) align the cutting plane with the surfaces of the 
apparent horizons in the region where the surfaces are 
closest to the cutting plane. We derive the control sys- 
tem in charge of setting the parameters <d 3 by the follow- 
ing condition: We demand that the angle between the 
mapped cutting plane and the x-axis at x l c be driven to 
the (weighted) average of the angles at which the mapped 
apparent horizons intersect the same x-axis. 

The input coordinates to the skew map are the coordi- 
nates in the distorted frame. For each horizon, we there- 
fore measure an angle in the distorted frame as follows: 
Let x 1 ^ be the distorted-frame ^-coordinate at which ap- 
parent horizon H intersects the segment connecting the 



2 The horizon centers are not invariant under the skew map 
because they do not necessarily lie on this line. This deviation, 
which is quantified in Eqs. (37)-(39), leads to a coupling of O(Q) 
with the translation, rotation, and scaling maps. 



centers of the excision boundaries. For j = 1, 2, we calcu- 
late the normal to the surface at this intersection point. 



This normal is projected into the x 



3-3 



plane. 



We define Q 3 H as the angle between the projected nor- 
mal and the x-axis. Thus, projecting the normal into the 
y = const, plane gives Q Z H and vice versa (see Fig. 6). We 
then compute a weighted average of the Q J H such that the 
horizon closer to the cutting plane has a larger effect on 
the skew angles, 



p)3 

^Avg 



w A J , 



Wb&r 



w A w{c A ) + w B w(c B y 

where wa,wb are averaging weights, 



wh = exp 



,Jnt 



(66) 



(67) 



Here C l H is the center of the excision boundary H . 

We measure O^vg m the distorted frame, and in this 
frame the cutting plane is always normal to the x-axis. 
Therefore, thinking about the desired result in the dis- 
torted frame, we see that the control system for the 
skew map should drive 0^ v „ to zero. This leads us to 
consider the following control error for the skew angles: 
Q 3 e = 0^ vg . Assuming that the function W is unity near 
the apparent horizons, we find that 



-1 + 0(0), 



(68) 



in agreement with Eq. (13). In deriving Eq. (68), it is 
helpful to observe that the partial derivative in Eq. (68) 
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is taken with the inertial-frame apparent horizon held 
fixed. For a fixed inertial-frame horizon, the only map 
that can change the shape (as opposed to merely the 
center) of the distorted-frame horizon (and thus 6a vg ) is 
the skew map. 

We do not use Q 3 e = 6 Avg for the entire evolution, 
however, because at early times, when the coordinate 
distance between the apparent horizons is larger than 
their combined radii, the skew map is not needed. Fur- 
thermore, the skew map can cause difficulties early dur- 
ing the run, especially during the "junk radiation" phase 
when the horizons are oscillating in shape. For this rea- 
son, we gradually turn on the skew map as the black 
holes approach each other. This is done by defining a 
roll-on function g that is zero when the horizons are far 
apart, and one when they are close together. This roll-on 
function is defined as 



1 - tanh 10 



„Int 



-Int 



°A ~~ U B 



(69) 



For values g < 10 -3 the skew map is turned off com- 
pletely; this is not strictly necessary, but it saves some 
computation. Given the function <?, we define the control 
error as 



Qi=50iv K -(i-9)0 i - 



(70) 



This control error drives the skew angles to zero when 
the black holes are far apart and drives 6 Avg t° zer0 as 
they approach each other. 



E. Shape control 



We define the shape map -Mshape as: 



that these discontinuities occur on subdomain bound- 
aries, they cause no difficulty with using spectral meth- 
ods. Around excision region B, /b^b, (?b, </>b) is chosen 
similarly. 

In previous implementations of the shape map [3], 
the functions /#/(?~ir, , </>#) were chosen to be smooth 
Gaussians centered around each excision boundary rather 
than to be piecewise linear functions. We find smooth 
Gaussians to be inferior for two reasons. The first is that 
piecewise linear functions are easier and faster to invert 
(the inverse map is required for interpolation to trial so- 
lutions during apparent horizon finding). The second 
is that for smooth Gaussians, it is necessary to choose 
the widths of the Gaussians sufficiently narrow so that 
the Gaussian for excision region A does not overlap the 
Gaussian for excision region B and vice versa, so that the 
maps and control systems for A and B remain decoupled. 
However, decreasing the width of the Gaussians increases 
the Jacobians of the map, producing coordinates that 
are stretched and squeezed nonuniformly. We found that 
this form of "grid-stretching" significantly increased the 
computational resources required to resolve the solution 
to a given level of accuracy. In other words, with smooth 
Gaussians we were forced to add computational resources 
just to resolve the large Jacobians. 

A map very similar to Eq. (71) is also described in [4]. 
The difference compared to this work is the choice of 
fH(T"H,6H,4>H), which corresponds to a different choice 
of domain decomposition. The control system used to 
choose the map parameters A|L(i) in [4] is also different 
than what is described here. 

We control the expansion coefficients \f m (t) of the 
shape map, Eq. (71), so that each excision boundary, as 
measured in the intermediate frame (t,x l ), has the same 
shape as the corresponding apparent horizon. In other 
words, we desire 



H 



fH(r H ,0H,<l>H) 



I'm 



(71) 

The index H goes over each of the two excised regions 
A and B, and the map is applied to the grid- frame 
coordinates. The polar coordinates {rH,^H,4>H) cen- 
tered about excised region H are defined by Eqs. (31)- 
(33), the quantities Yg m (B h , <j> jj) are spherical harmon- 
ics, and Xf m (t) are expansion coefficients that param- 
eterize the map near excision region H. The function 
fii{rH,&H,<l>H) is chosen to be unity near excision re- 
gion H and zero near the other excision region, so that 
the distortion maps for the two black holes are decou- 
pled. Specifically, </>#) is determined as il- 
lustrated in Fig. 4. For excision region A in the fig- 
ure, Ja^a? 4> a) = 1 everywhere inside the magenta 
sphere, it is zero everywhere outside the surface shown 
in red, and it falls linearly to zero between these two sur- 
faces. This means that the gradient of Ja^a, &a, 4>a) is 
discontinuous on the red surface, the magenta surface, 
and the blue surfaces in the figure. Because we ensure 



^AH 

(^ah) 



(H3b) 



(72) 



where for brevity we have dropped the index H that la- 
bels the excision boundary. Here the angle brackets mean 
averaging over angles, tah(0, 4>) is the radial coordinate 
of the apparent horizon defined in Eq. (35), and teb ifi, 4>) 
is the radial coordinate of the excision boundary, which 
from Eq. (71) can be written 



^eb = r EB - ^2 Y lm (6, 4>)X im (t). 



(73) 



Here teb is the radius of the (spherical) excision bound- 
ary in the grid frame. In deriving Eq. (73) we have used 
the relations A^cutx = 1, f(r,0,<j)) = 1, 6 = 6, and 
</> = </>, which hold on the excision boundary. 
Combining Eqs. (35), (72), and (73) yields 



teb - Y 0Q X Q0 (t) 



SqoYqo 



(74) 
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which we can satisfy by demanding that 

Mm{t) + b em g =0, I > 0. (75) 

Therefore, given an apparent horizon and given a value 
of Aoo, a control system can be set up for each (I, m) pair 
with £ > 0, and the corresponding control errors are 

n \ m c r EB -Y 00 X 00 (t) 

Qtm = -Mm{t) - btm g~ ~ , « > 0. (76) 

Driving these Qi m to zero produces an excision boundary 
that matches the shape of the corresponding apparent 
horizon. 

Eq. (76) determines Xi m {t) only for I > 0, and 
Eqs. (72)-(76) can be satisfied for arbitrary values of the 
remaining undetermined map coefficient Xoo(t). Determi- 
nation of Xoo(t) is complicated enough that it is described 
in its own section, Sec. V. 

Note that Eq. (76) does not satisfy Eq. (14) because 
dQim/ dXoo — Sim/Soo, which does not vanish even if 
all the control errors are zero. For small distortions, this 
coupling between Aoo and Qg m seems to cause little dif- 
ficulty However, at times when the shapes and sizes of 
the horizons change rapidly (e.g., during the initial "junk 
radiation" phase, after the transition to a single excised 
region when a common apparent horizon forms, and after 
rapid gauge changes), simulations using Eq. (76) exhibit 
relatively large and high-frequency oscillations in Qi m 
that usually damp away but occasionally destabilize the 
evolution. Construction of a control system in which all 
Qim are fully decoupled will be addressed in a future 
work. 



F. CutX 

The map Mcutx applies a translation along the grid- 
frame s-axis in the vicinity of the black holes, but with- 
out moving the excision boundaries (or the surrounding 
spherical shells) themselves. The action of the map is 
shown in Fig. 7. 

The goal of this map is to allow for a slight motion 
of the cutting plane towards the smaller excision bound- 
ary. This is important for binary black hole systems with 
mass ratio q > 8. As such a binary gets closer to merger, 
the inertial-frame coordinate distance between each ex- 
cision boundary and the cutting plane decreases. Even- 
tually, this distance falls to zero for the larger excision 
boundary, producing a coordinate singularity in which 
the Jacobian of one of the other maps (often the shape 
map) becomes infinite. By pushing the cutting plane to- 
wards the smaller excision boundary, the map A^cutx 
avoids this singularity. Even for evolutions in which the 
inertial-coordinate distance between the larger excision 
boundary and the cutting plane remains finite but be- 
comes small, the map A^cutx prevents large Jacobians 
from developing and thus increases numerical accuracy 




FIG. 7: Two-dimensional projection of the domain decom- 
position near the two black holes A and B. The function p 
in Mcutx is zero on the magenta boundaries, one on the red 
boundaries, and goes linearly between zero and one in the "ra- 
dial" direction between these boundaries. The gradient of p is 
discontinuous across the solid magenta, red, and blue bound- 
aries. Dotted lines show the subdomain boundaries under the 
action of A^cutx- 

(because it is no longer necessary to add computational 
resources to resolve the large Jacobians). 

Figure 8 shows the domain decomposition in the iner- 
tial frame for a binary black hole simulation with q = 8. 
The top panel shows the case without A^cutx, and the 
compressed grid near the larger excision boundary is evi- 
dent. The lower panel shows the case with A^cutx, which 
removes the extreme grid compression. 

The map .Mcutx is written as 

x° ^ x° + p(x l )F(t), (77) 
x j ^x j , j = 1,2, (78) 

where F(t) is adjusted by a control system. 

We choose p(x z ) to be zero within either of the spheri- 
cal shell regions, i.e., inside the magenta spheres around 
A and B in Fig. 7, and it is also zero outside the outer 
magenta sphere in the same figure. We set p(x l ) = 1 on 
the solid red boundaries in Fig. 7, and everywhere be- 
tween the two solid red boundaries that enclose excision 
boundary B. Every other region on the grid is bounded 
by a smooth red boundary on one side and a smooth ma- 
genta boundary on the other; in these regions p varies 
linearly between zero and one. The solid blue bound- 
aries are locations (in addition to the solid red and solid 
magenta boundaries) in which the gradient of p is dis- 
continuous. Full details for the calculation of p can be 
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FIG. 8: Snapshot of the grid, viewed in the inertial frame, 
from a binary black hole evolution with Ma = 8Mb- Top: 
without Alcutx- Bottom: with A^cutx- The gray areas are 
excised regions. In the top panel, the grid is compressed as 
the larger excision boundary approaches the cutting plane; 
this is especially evident in the long narrow subdomains im- 
mediately adjacent to the larger excision boundary. Using 
-Mcutx relieves the grid compression. 



found in Appendix A. 

As with the skew map, the map A^cutx is inactive for 
most of the inspiral. Let xff c be the distorted-frame x- 
coordinate of the intersection of the excision boundary 
H with the segment connecting the centers of the exci- 
sion boundaries. This is similar to x 1 ^ defined earlier; 



Cj-f* refers to a point on the appar- 



the difference is that x\ 

j,\x re f ers t a point on the excision 



ent horizon and x" H 

boundary. The quantity x^ c is time-dependent because 
it depends on the shape map. 

The map A^cutx is turned off completely as long as 

„0 ^.Exc 



^Exc 



-c% 



< 



1 



(79) 



H 



where C l H are the excision boundary centers, and x % c are 
the coordinates of the intersection of the cutting plane 
and the segment connecting the centers of the excision 
boundaries, as introduced in Sec. IV D. 

Let to be the coordinate time at which the map A^cutx 
is activated. We designate the target position x = xt of 
the cutting plane as 

Ax a • y Exc 1 A - - ~- Exc 
xt = - 



Ax A 



Ax B 
Ax B 



(80) 



where 



Ax H 



l<r Exc 
\ X H 



r*o I 



Recall that ir Exc is time-dependent (as it is measured 
in the distorted frame), so we save the value of xt as 
calculated at the time when the map .Mcutx is activated, 
rather than recalculating xt at every measurement time. 

Now that we have a target xt, we could designate 
xt — x*q as the target value to the function F(t) by set- 
ting Qp = —F + xt~ Xq and have the control system 
drive Qf to zero, thus driving the x-coordinate of the 
cutting plane to xt- However, because we turn on the 
CutX control system suddenly at time to, we must be 
more careful. Turning on any control system suddenly 
will produce some transient oscillations, unless the con- 
trol error and its relevant derivatives and integrals are 
all initially zero. In the case of the CutX map, which 
is turned on during a very dynamic part of the simula- 
tion when excision boundaries need to be controlled very 
tightly, these oscillations can prematurely terminate the 
run by, for example, pushing an excision boundary out- 
side its accompanying apparent horizon. 

To turn on A^cutx gradually, we replace xt by a new 
time-dependent target function T(t) that gradually ap- 
proaches Xt at late times but produces a control error Qp 
with Qp = dQpjdt = at the activation time t = to- 
We start by estimating the time tfp at which x^f c will 
reach the cutting plane, where H is either A or B. This is 
done by the method described in Appendix B. We then 
let t be the minimum of and t^ c . The smooth 
target function T(t) is then defined by 



T(t;t ,%T,t 



XC\ 



xt + exp 



x° c - x T + F(t ) + (t-io) 



(t - tp) 
0.3 • t xc 

dF(t) 



at 



t=t 



Designating T the target for F(t) 



Q F = T-F-x° c 



leads to 



Qp\t=t — 0, 



dQi 



dt 



(82) 



(83) 



(84) 



t=t 



so at the activation time the control system does not 
produce transients. Furthermore, in the limit of small 
Qf, 



x° c + F(t)\ t=t XC w X T , 



(85) 



H = A,B. 



i.e., by the time the excision boundary would have 
touched the cutting plane (and formed a grid singularity), 
the smooth target function T(t) has approached Xt, and 
the cutting plane will have reached its designated target 
location, xt- As the run proceeds, the behavior of the 
cutting plane is determined by the map A^cutx while the 
motion of the excision boundaries is determined by gauge 
dynamics and the behavior of the other control systems. 
We continue to monitor the distance between the cutting 
plane and the excision boundaries, and if it is predicted 
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to touch within a time less than Td/0.15, the A^cutx con- 
trol system is reset, constructing a new target function 
T(t;t ,XT,t xc ), where to is the time of the reset, and 
xx,t xc are also recalculated based on the state of the 
grid at this reset time. 

Each time a new target function T is constructed, the 
damping time of the A^cutx control system is set to 
bet XG /2. 

The control systems responsible for A^cutx and 
•Mshapo are decoupled, as Mcutx controls the location 
of the the cutting plane, leaving the excision boundary 
unchanged, while A^shape controls the shape of the exci- 
sion boundary, leaving the location of the cutting plane 
unchanged. Recall that A^cutx and A'lshapc define the 
mapping from the grid frame to the distorted frame. As 
the apparent horizons are found in the distorted frame, 
and the other maps only depend upon measurements of 
the horizons, A^cutx and A^shapc are decoupled from the 
other maps. 



V. SIZE CONTROL 

In this section we discuss how we control the spherical 
part of the map given by Eq. (71), namely, the coefficients 
Aqo for each excision boundary H. We apply the same 
method to each excision boundary, so for clarity, in this 
section we again drop the index H from the coefficients 
Xoo and Sf m , and from the coordinates (th , Oh-, <Az? )• 



A. Direct size control 

One possibility is to demand that the excision bound- 
ary radius is some fraction of the apparent horizon radius, 
i.e., 



r E B = crr A H 
where a < 1. This leads to the control error 



Q 



teb — vSqoYoo 



1™ 



A 



no ■ 



(86) 



(87) 



A more sophisticated version of this control system was 
employed in [4]. 

Although this control system is effective in keeping 
the excision boundary tied to the apparent horizon, it 
is not very robust in binary simulations, especially when 
the black holes are rapidly changing in size (because of 
physics or gauge). This is because black hole excision 
requires conditions on the characteristic speeds of the 
system, and if these conditions are not enforced they are 
likely to be violated. 



B. Characteristic speed control 

The minimum characteristic speed at each excision 
boundary is given by 



dx 1 



where a is the lapse, /3 l is the shift, and hi is the normal 
to the excision boundary pointing out of the computa- 
tional domain, i.e., towards the center of the hole. Here 
(i, x l ) are the inertial frame coordinates. The first two 
terms in Eq. (88) describe the coordinate speed of the 
ingoing (i.e., directed opposite to hi) light cone in the in- 
ertial frame, and the last term accounts for the motion of 
the excision boundary (which is fixed in the grid frame) 
with respect to the inertial frame. 

In our simulations, we impose no boundary condition 
whatsoever at each excision boundary. Therefore, well- 
posedness requires that all of the characteristic speeds, 
and in particular the minimum speed v, must be non- 
negative; in other words, characteristics must flow into 
the hole. In practice, we indeed find that if v becomes 
negative, the simulation crashes soon afterward. This 
can occur even when using the control system described 
by Eq. (87) and variations thereof. 

Therefore, we would like to control Aoo in such a way 
that v remains positive. We start by writing v in a way 
that separates terms that explicitly depend on Aoo from 
terms that do not. To do this we expand the derivative 
in the last term of Eq. (88) as 

dx 1 dx 1 dx a 
~dt ~ d~i°-~dt 

dx 1 di dx 1 dxi 
~ Itidt + d&lh 

dx 1 dx 1 dxi 
~ ~dT + d&~dt' 



(89) 



where a in the first line of Eq. (89) is a four-dimensional 
spacetime index, and the last line of Eq. (89) follows from 
dt/dt — 1. Inserting this into Eq. (88) yields 



dx 1 



dx 1 



(90) 



where x l is the frame which is obtained from the grid 
frame by applying the distortion map; see Eq. (29). 
We then use Eq. (71) to rewrite this as 



-a — f3 l ni 



dx 1 



tra 



(91) 



where we have used the relations f(r, 0,(f>) — 1 and 
A^cutx = 1 when evaluating the distortion map 
■M Distortion on the excision boundary. 
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By combining all the terms that do not explicitly de- 
pend on Ago into a quantity vq, we obtain 



v = v + hi—YonXoo- 
r 



(92) 



Thus, the characteristic speed v can be thought of as con- 
sisting of two parts: one part, Uq, that depends on the 
position and shape of the excision boundary and the val- 
ues of the metric quantities there, and another part that 
depends on the average speed of the excision boundary 
in the direction of the boundary normal. 

We now construct a control system that drives the 
characteristic speed v to some target speed vt- This is a 
control system that controls the derivative quantity Aoo, 
as opposed to directly controlling the map quantity Aoo- 
We choose 



Q = (min(w) - v T )/(-E), 



where 



S — hi — Yoo, 
r 



(93) 



(94) 



the minimum is over the excision boundary, and the angle 
brackets in Eq. (93) refer to an average over the excision 
boundary. Note that S < because hi points radially 
inward and x l /r points radially outward; this means that 
Q = — A o, in accordance with our normalization choice 
for a control system on Aoo- 

As in our other controlled map parameters, we demand 
that Aoo is a function with a piecewise-constant second 
derivative. It is then easy to construct Aoo as a function 
with a piecewise-constant third derivative. 

The control system given by Eq. (93) with a hand- 
chosen value of vt has been used successfully [24, 25] in 
simulations of high-spin binaries. Fig. 9 illustrates why 
characteristic speed control is crucial for the success of 
these simulations. 



C. Apparent horizon tracking 

The characteristic speed control described above has 
the disadvantage that it requires a user-specified target 
value vt- If vt is chosen to be too small, then small 
fluctuations (due to shape control, the horizon finder, or 
simply numerical truncation error) can cause the charac- 
teristic speed to become negative and spoil the simula- 
tion. 

If vt is chosen to be too large, the simulation can also 
fail. To understand why, recall that characteristic speed 
control achieves the target characteristic speed by mov- 
ing the excision boundary and thus changing the velocity 
term in Eq. (88). So if v > vt the control system moves 
the excision boundary radially inward, and if v < vt the 
control system moves the excision boundary radially out- 
ward. If vt is too large, the control system can push the 
excision boundary outward until it crosses the apparent 
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FIG. 9: Minimum characteristic speed on one of the excision 
surfaces of an binary with equal aligned spins of \ — 0.97, 
shown just before merger [25]. Two cases are shown: one 
without characteristic speed control (solid) and one that is 
restarted from the uncontrolled run at t = 6255Af with char- 
acteristic speed control turned on and vt = 5 x 10 (dashed). 
The uncontrolled run crashes because of negative character- 
istic speeds (i.e., ill-posedness) shortly after t — 6290M, but 
the controlled run continues through merger and ringdown. 



horizon. This halts the evolution because the apparent 
horizon can no longer be found. 

One way to prevent the excision boundary from cross- 
ing the horizon is to drive the excision boundary to some 
constant fraction of the horizon radius, or in other words, 
drive the quantity d/dt{Ar) to zero, where 



Ar = 1 



(^ah) 



(95) 



is the relative difference between the average radius of 
the apparent horizon (in the intermediate frame) and the 
average radius of the excision boundary. Using Eqs. (73) 
and (35), we can write 



d A 

— At 
dt 



Aoo 
Sao 



S, 



(1 - Ar) . 



(96) 



on 



and therefore a control system that adjusts Aoo to achieve 
d/dt(Ar) = can be obtained by defining 



S 00 (Ar-l)-A 



(97) 



A slight generalization of this control system can be 
obtained by demanding that d/dt(Ar) — rariftj where 
^drift is some chosen constant, 



#oo (Ar - 1) - Aoo + Soo?\Mft- 



(98) 



We have not found the control systems defined by 
Eqs. (97) and (98) to be especially useful on their own. 
One drawback of these systems is that they do not pre- 
vent the minimum characteristic speed v from becoming 
negative. Instead, we use the horizon tracking control 
systems described here as part of a more sophisticated 
control system discussed in the next section. 
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D. Adaptive switching of size control 

Here we introduce a means of controlling Aoo that com- 
bines the best features of characteristic speed control and 
horizon tracking. The idea is to continuously monitor the 
state of the system and switch between different control 
systems as the evolution proceeds. 

At fixed intervals during the simulation which we call 
"measurement times," we monitor the minimum char- 
acteristic speed v and the relative distance between the 
horizon and the excision boundary Ar. The goal of the 
control system is to ensure that both of these values re- 
main positive. 

Consider first the characteristic speed v. Using the 
method described in Appendix B, we determine whether 
v is in danger of becoming negative in the near future, 
and if so, we estimate the timescale t v on which this will 
occur. Similarly, we also determine whether Ar will soon 
become negative, and if so we estimate a corresponding 
timescale TA r - Having estimated both t v and TAr, we 
use these quantities to determine how to control Aoo- In 
particular, we switch between horizon tracking, Eq. (97), 
and characteristic speed control, Eq. (93), based on t v 
and TAr- 

We do this by the following algorithm, which favors 
horizon tracking over characteristic speed control unless 
the latter is essential. Assume that the control system for 
Aoo is currently tracking the horizon, Eq. (97), and that 
the current damping timescale is some value T d . If v is in 
danger of crossing zero according to the above estimate, 
and if t v < and t v < ta t , we then switch to charac- 
teristic speed control, Eq. (93), we set vt = 1.01k where 
v is the current value of the characteristic speed (the 
factor 1.01 prevents the algorithm from switching back 
from characteristic speed control to horizon tracking on 
the very next time step) , and we reset the damping time 
Td equal to t v . Otherwise we continue to use Eq. (97), 
resetting r d = TAr if i~Ar < T d . 

Now assume that the control system for Aoo is con- 
trolling the characteristic speed, Eq. (93). If v is not in 
danger of crossing zero, so that we no longer need ac- 
tive control of the characteristic speed, then we switch to 
horizon tracking, Eq. (97), without changing the damp- 
ing time T d . If v is in danger of crossing zero, and if Ar is 
either in no danger of crossing zero or if it will cross zero 
sufficiently far in the future such that TAr > &\Td, then 
we continue to use Eq. (93) with Td reset to min(rd, t v ) so 
as to maintain control of the characteristic speed. Here 
(7 1 is a constant that we typically choose to be about 20. 

The more complicated case occurs when both v and 
Ar are in danger of crossing zero and t Ar < &iT~d- In 
this case, we have two possible methods by which we at- 
tempt to prevent Ar from becoming negative. We choose 
between these two methods based on the values of ta v 
and Td, and also by the value of a quantity that we call 



the comoving characteristic speed: 

ai- - d ^ 

v c = -a- p m - rii—^ 
at 

+ n t —l r o^oo(Ar - 1) + V Y im (6, <t>)\ im (t) ) . 

(99) 

This quantity is the value that the characteristic speed v 
would have if horizon tracking (Eq. (97)) were in effect 
and working perfectly. Eq. (99) is derived by assuming 
Q = in Eq. (97), solving for Aoo, and substituting this 
value into Eq. (91). The reason to consider v c is that 
for min(i! c ) < (which can happen temporarily during a 
simulation), horizon tracking is to be avoided, because 
horizon tracking will drive min(w) towards a negative 
value, namely min(i; c ). 

The first method of preventing Ar from becoming neg- 
ative is to switch to horizon tracking, Eq. (97); we do 
this if min(i! c ) > and if t Ar < C2i~d, where 02 is a 
constant that we typically choose to be 5. The second 
method, which we employ if (J2 T d < r Ar < &i T d or if 
min(v c ) < 0, is to continue to use characteristic speed 
control, Eq. (93), but to reduce the target characteris- 
tic speed vt by multiplying its current value by some 
fraction rj (typically 0.25). By reducing the target speed 
vt, we reduce the outward speed of the excision bound- 
ary, and this increases TAr- The reason we use the lat- 
ter method for min(u c ) > is that the former method, 
horizon tracking, will drive v towards v c and we wish to 
keep v positive. The reason we use the latter method for 
&2Td < TAr < &ii~ d is that in this range of timescales the 
control system is already working hard to keep v posi- 
tive, so therefore a switch to horizon tracking is usually 
followed by a switch back to characteristic speed con- 
trol in only a few time steps, and rapid switches between 
different types of control make it more difficult for the 
control system to remain locked. The constants r/, a±, 
and <72 are adjustable, and although the details of the 
algorithm (i.e., which particular quantity is controlled 
at which particular time) are sensitive to the choices of 
these constants, the overall behavior of the algorithm is 
not, provided 1 < 02 < <xi. 

The overall behavior of this algorithm is to cause Aoo 
to settle to a state in which both v > and Ar > for a 
long stretch of the simulation. Occasionally, when either 
v or Ar threatens to cross zero as a result of evolution 
of the metric or gauge, the algorithm switches between 
characteristic speed control and horizon tracking several 
times and then settles into a new near-equilibrium state 
in which both v > and Ar > 0. 

The actual value of Ar and of v in the near-equilibrium 
state is unknown a priori, and in particular this algorithm 
can in principle settle to values of Ar that are either 
small enough that control system timescales must be very 
short in order to prevent Ar from becoming negative, or 
large enough that excessive computational resources are 
needed to resolve the metric quantities deep inside the 
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horizon. To prevent either of these cases from occur- 
ring, we introduce a constant correction velocity f corr , 
a nominal target value Ar ta rget, and a nominal target 
characteristic speed ^target- We currently choose these 
to be f corr = 0.005, Artarget = 0.08, and target = 0.08. 
These quantities are used only to decide whether to re- 
place Eq. (97) by Eq. (98) when the algorithm chooses 
horizon tracking, as described below; these quantities are 
not used when the algorithm chooses characteristic speed 
control. 

First we consider the case where Ar is too large during 
horizon tracking. If Ar > l.lAr targct , then we replace 
Eq. (97) with Eq. (98), and we use Thrift = —r C on- We 
continue to use r drift = — ^corr until Ar < Ar ta rgot, at 
which point we set f drift = 0. 

The case where Ar is too small during horizon track- 
ing is more complicated because we want to be care- 
ful so as to not make the characteristic speeds decrease, 
since it is more difficult to recover from decreasing char- 
acteristic speeds when Ar is small. In this case, if 
Ar < 0.9Ar tar got, v < 0.9v t argct, and (n k d k v c ) > 0, 
then we replace Eq. (97) with Eq. (98), and we use 
f'drift = w corr . Here v c is the comoving characteristic 
speed defined in Eq. (99) and the angle brackets refer 
to an average over the excision boundary. The condition 
on the gradient of v c ensures that v c will increase if Ar 
is increased; otherwise a positive rdrift will threaten the 
positivity of the characteristic speeds. We continue to use 

Adrift — ^corr 

until either v > 0.9z;target, Ar > 0.9Ar ta rget, 
(n k dkV c ) < 0, or if v will cross ^target or Ar will cross 
Artarget within a time less than the current damping 
time Td (where crossing times are determined by the pro- 
cedure of Appendix B). When any of these conditions 
occur, we switch to rdrift = 0, and we prohibit switch- 
ing back to a positive rdrift until either v > 0.99utargct or 
Ar > 0.99 Ar t argct j or until both v and Ar stop increas- 
ing in time. These last conditions prevent the algorithm 
from oscillating rapidly between positive and zero values 
of r dr ift- 



VI. CONTROL SYSTEMS FOR EFFICIENCY 

Although the most important use of control systems 
in SpEC is to adjust parameters of the mapping between 
the inertial and grid coordinates, another situation for 
which control systems are helpful is the approximation of 
functions that vary slowly in time, are needed frequently 
during the simulation, but are expensive to compute. 

For example, the average radius of the apparent hori- 
zon, rAH, is used in the control system for Aoo, which is 
evaluated at every time step. Evaluation at each time 
step is necessary to allow the control system for Aoo to 
respond rapidly to sudden changes in the characteristic 
speed or the size of the excision boundary, as may occur 
after regridding, after mesh refinement changes, or when 
other control systems (such as size control) are temporar- 
ily out of lock. 



However, computing the apparent horizon (and thus its 
average radius rAH) at every time step is prohibitively 
expensive. To significantly reduce the expense, we de- 
fine a function with piecewise-constant second derivative, 
r AH X (*)' we define a control error 

Q = r Aii - rfi\ (100) 

and we define a control system that drives Q to zero. We 
then pass r A jj X (i) instead of rAH to those functions that 
require the average apparent horizon radius at each time 
step. The control error Q, and thus the expensive com- 
putation of the apparent horizon, needs to be evaluated 
only infrequently, i.e., on the timescale on which the av- 
erage horizon radius is changing, which may be tens or 
hundreds of time steps. 



VII. SUMMARY 

In simulations of binary black holes, we use a set 
of time-dependent coordinate mappings to connect the 
asymptotically inertial frame (in which the black holes 
inspiral about one another, merge, and finally ringdown) 
to the grid frame in which the excision surfaces are sta- 
tionary and spherical. The maps are described by pa- 
rameters that are adjusted by a control system to follow 
the motion of the black holes, to keep the excision surface 
just inside the apparent horizons of the holes, and also to 
prevent grid compression. We take care to decouple the 
control systems and choose stable control timescales. 

The scaling, rotation, and translation maps are used 
to track the overall motion of the black holes in the in- 
ertial frame. These are the most important maps during 
the inspiral phase of the evolution, as the shapes of the 
horizons remain fairly constant after the relaxation of the 
initial data. As the binary approaches merger, the hori- 
zon shapes begin to distort, and shape and size control 
start to become important. Shape and size control are 
especially crucial for unequal-mass binaries and for black 
holes with near-extremal spins. In the latter case, the 
excision surface (which must remain an outflow bound- 
ary with respect to the characteristics of the evolution 
system) needs to be quite close to the horizon. 

In order to decouple the shape maps of the individual 
black holes, our grid is split by a cutting plane at which 
the shape maps reduce to the identity map. As the black 
holes merge, a skew map is introduced in order to align 
the cutting plane with the excision surfaces (see Fig. 5) 
and to minimize grid compression between the two black 
holes. 

Finally, the implementation of CutX control was nec- 
essary to complete the merger of high mass-ratio con- 
figurations. In these systems, the excision boundaries 
approach the cutting plane asymmetrically as the black 
holes merge, thus compressing the grid between the cut- 
ting plane and the nearest excision boundary. CutX con- 
trol translates the cutting plane to keep it centered be- 



18 



tween the two excision boundaries, alleviating this grid 
compression. 

We find that we need all of the maps described in this 
paper with a sufficiently tight control system in order to 
robustly simulate a wide range of the parameter space of 
binary black hole systems [4, 5, 24]. Many of these maps 
and control systems are also used in our simulations of 
black hole - neutron star binaries [26-28]. 



Appendix A: Details on computing the CutX map 
weight function 

As discussed in Sec. IV F, we use the map .Mcutx 
to control the location of the cutting plane in the last 
phase of the merger, as the distance between the exci- 
sion boundaries and the cutting plane becomes small. 
The value of the weight-function p(x l ) in Eq. (77) is zero 
in the spherical shells describing the wave zone and in the 
shells around the excision boundaries (see Fig. 7). The 
weight function is unity on the cutting plane, and on the 
"cut-sphere" surfaces around either excision boundary. 
For the smaller excision surface we have two such cut- 
sphere surfaces and the value of the weight function is 
one within the region enclosed by these two cut-spheres 
and the cutting plane. The weight function transitions 
linearly from zero to one (from the magenta curves to the 
red curves on Fig. 7). In these transitional regions the 
value of p{x l ) is computed as described below. 

Consider the region spanning the volume between the 
spherical shells around the excision boundary H (with 



H = A,B) and the cutting plane x 



(This is 



referred to in the code as the M region). In order to 
calculate the value of p(x 1 ) in this region, we shoot a 
ray from C l H in the direction of x' . This ray intersects 
both the outer spherical boundary of the shells and the 
cutting plane. The ir-coordinate of the intersection point 
with the shells will be 



(x° - C%)R H 



c 



(Al) 



H) 



where Rr is the outer radius of the spherical region 
around C l H . Then, for a point x % in the M region p 



is given by 



p = (x° - xm)/(x° c - x M )- 



(A2) 



All other regions affected by the map Mcutx are de- 
limited on both ends by spheres. Here we proceed as 
follows: Let x l P denote the "center of projection," i.e., 
the point toward which the unmapped radial grid lines 
of the given zone converge, and let x % s denote the center 
of a sphere of radius R. We will need to solve for the 
parameter ?7(a;s, £p, x, R) defined by 



(xp + T](x l - x P ) - x x ^\ 



= R 2 



(A3) 



We want to use the positive solution of the quadratic 
system, 



n 



B + VB 2 - C, 



where 



B = 



C 
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(A4) 

(A5) 
(A6) 



Given two spheres delimiting the region where we need 
to compute p, we proceed as follows. On one of these 
spheres we want p — so we solve for Eq. ( A4) , obtaining 
t]q. We repeat this procedure to obtain r]\ on the spherical 
surface where we want p = 1. These values of 77 define 
intersection points of the two spherical surfaces and the 
line connecting x l with the projection center, x P . We 
denote the distance between these points and x P by rg 



and 7"i , while the distance between x % 
p is given by 



and 



p = (r- r )/(ri - r ). 



is r. Then 



(A7) 



On either side of the cutting plane, computing p between 
the inner shells and cut-sphere is done using C l H as the 
projection center. The value of p vanishes on the inner 
sphere and it equals one on the outer, cut-sphere. Next, 
when computing the value of p between the cut-sphere 
and the larger, outer sphere, the projection center is x l c . 
The value of p is one on the cut-sphere and it vanishes 
on the outer, larger sphere. 

The input coordinates to VWcutx are the coordinates 
x l that are obtained from the grid coordinates x l by the 
shape map, x 1 — Ms\ia,pc{ x% ) ■ As seen above, computa- 
tion of the weight function at any point requires knowl- 
edge of the region containing that point, which in turn re- 
quires knowledge of its grid- frame coordinates x l . There- 
fore, to compute the weight function one must first com- 



pute x 



M 



Shape 



x l ). When computing the forward 



CutX map, this inverse is done only once. However, the 
situation is more complicated when computing the in- 
verse CutX map, because the grid-frame coordinates de- 
pend upon the inverse CutX map itself, which depends 
on the weight function. This leads to an iterative algo- 
rithm where we must call the inverse map function of 
■M shape from each iteration of the inverse map function 
of A^cutx- This could easily become an efficiency bot- 
tleneck. We mitigate this by keeping Mcutx inactive 
for most of the run, only activating it when it becomes 
crucial in order to avoid a singular grid map. 



Appendix B: Estimation of crossing times 

Several of the control systems described here are de- 
signed to ensure that some measured quantity remains 
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positive. For example, in Sec. VD, we demand that both 
the characteristic speed v at the excision boundary and 
the difference Ar between the horizon radius and the 
excision boundary radius remain positive. Similarly, in 
Sec. IV F we demand that each excision boundary does 
not cross the cutting plane. 

Part of the algorithm for ensuring that some quantity 
q remains positive is estimating whether q is in danger 
of becoming negative in the near future, and if so, esti- 
mating the timescale r on which this will occur. In this 
Appendix we describe a method of obtaining this esti- 
mate. 

We assume the quantity q is measured at a set of (not 
necessarily equally-spaced) measurement times, and we 
remember the values of q at several (typically 4) previous 
measurement times. At each measurement time to we fit 
these remembered values to a line q(t) — a + b(t — to). 
The fit gives us not only the slope b and the intercept 
a, but also error bars for these two quantities da and 
8b. Assuming that the true q(t) lies within the error 
bars, the earliest time at which q will cross zero is t = 
t Q + (—a+5a)/(b—5b), and the latest time is t = t + (—a— 
5a) jib + Sb). If both of these times are finite and in the 
future, then we regard q as being in danger of crossing 
zero, and we estimate the timescale on which this will 



happen as r = —a/b. 
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